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A COMPARISON OF IMPLICIT NUMERICAL 


METHODS FOR SOLVING THE TRANSIENT 

SPHERICAL DIFFUSION EQUATION 

By Donald M. Curry 
Lyndon B. Johnson Space Center 


SUMMARY 


Comparative numerical temperature results obtained by using two implicit 
finite-difference procedures for the solution of the transient spherical heat- 
conduction equation are presented. The strongly implicit procedure is com- 
pared to the more standard alternating-direction implicit procedure by using 
a two-dimensional solid spherical model. The numerically generated tempera- 
ture results obtained by using the strongly implicit procedure and the 
alternating-direction implicit procedure » re compared with exact solutions to 
assess the relative accuracy and efficiency of the two numerical methods. 
Special attention was given to the solution in the regions of singularities 
associated with the governing partial differential equation. For the examples 
solved, the numerical results obtained by a modified version of the strongly 
implicit procedure and by the alternating-direction implicit procedure are 
in close agreement with the exact solution. 


INTEC-LiUCTION 


Numerous authors have discussed the various numerical methods available 
for solving the transient diffusion equation. Solutions to the diffusion 
equation by means of numerical methods are required for a wide variety of 
design/development problems associated with the aerospace, petroleum, and 
chemical industries. 

Trent and Welty (ref. l) presented a good summary of numerical methods 
for solving transient-heat-conduction problems. However, they did not include 
discussion of a recently developed iterative technique (Stone, ref. 2) 
called the strongly implicit procedure (SIP). The SIP was shown to have 
several advantages over other implicit numerical techniques in solving large 
sets of algebraic equations that arise in the approximate solution of multi- 
dimensional partial differential equations. Weinstein et al. (ref. 3) have 
used the SIP successfully to solve systems of equations arising in multi- 
phase, two-dimensional reservoir flow problems. The SIP has been used by 
Curry (ref. L) in the solution of two-dimensional heat and mass transfer in 
porous media. Steen and All (ref. 5) compared the SIP algorithm with the 


more conventional implicit method in the solution of the nonlinear partial 
differential equation for the flow of a real gas in two dimensions. However, 
few two-dimensional numerical solutions of the transient-heat-conduction 
equation for both spherical and cylindrical coordinates can be found in the 
literature. Albasiny (ref. 6) presented an implicit numerical solution for 
a cylindrical heat-conduction problem, including the effects of the singularity 
at the center of the solid. Kee (ref. 7) developed a finite-difference algo- 
rithm for the diffusion equation for a solid sphere. 

In this report, the SIP is compared with the more conventional 
alternating-direction implicit procedure (ADIP) (ref. 8) by us'ng a two- 
dimensional spherical heat-conduction model. The temperature results obtained 
are compared to exact solutions of the spherical heat-conduction equation for 
various boundary conditions. Attention is given to the adequacy of the finite- 
difference representation in the neighborhood of the singularities located at 
the geometrical center, r * 0, and along the boundaries, 4> * 0, n. 


SYMBOLS 

A, B, C, D, E, Q parameters known from previous time level and previous 
iteration 

C specific heat 

P 

k thermal conductivity 

m, n iterative variables used in equation (lU) 

q* * * volumetric heat source (sink) 

R outer sphere radius 

\ 

T temperature 

T' unknown temperature in difference equations 

t time 

r, $, 0 spherical space coordinates 

x, y rectangular space coordinates 

Y iteration parameter 

p density 


2 


Subscripts 


1 

J 

x 

y 

Superscript 


♦-direction node location 
r-direction node location 
x-di rection 
y-direction 

indicates parameter at .ext time step or iteration 


THEORETICAL MODEL 


The transient-heat-conduction equation in spherical coordinates, with the 
assumption of constant thermophysical properties, is given as 


pc 


3T . 1 

at ■ k [r 


1 a 2 (rT) 

, 1 3 

/,). . 3T\ 

2 " 

1 3T 

»/ 

2 . . at 

r sin ♦ 

V 1 * *♦/ 

2 2 

r sin ♦ 30 ^ 


(1) 


where p is density; r, ♦, and 0 are spherical space coordinates, defined 
in figure 1; T is temperature; t is time; k is thermal conductivity; C 
is specific heat; and q ' * * is volumetric heat source. 

If the temperature field has azimuthal symmetry, then 



(2) 


Equation (l) can then be written in two dimensions as 



3T 

at 


a T 2k 3T 

V r 9r 


r^sin ♦ 


3 ? I 81 " ♦ fi) 


( 3 ) 


This two-dimensional unsteady heat conduction in the spherical domain is 
bounded by 


with the boundary and initial conditions as 


T(R, ♦, t) - fjU, t) 


(«•) 


where R is the outer sphere radius. 

In formulating the boundary conditions, it should be noted that equation 
(3) is singular at r ■ 0 and for $ ■ 0, *. The boundary condition 
represented by equation ( U ) permits a sphere with a variable surface tempera- 
ture from ♦ ■ 0 to $ * a (i.e., a sphere that is hot at the top and cold 
at the bottom). This variation obviously will result in a temperature 
gradient at r * 0. For this analysis, it is assumed that 


3T 

3r 


0 , 


r 


0 


(5) 


Equation (5) is strictly true only on $ * n/2. 1 

On the assumption of symmetry along $ ■ 0, n, then 

~ ■ 0, * ■ 0, n (6) 

The initial condition is 

T(r, ♦, 0) ■ f 2 (r, ♦) (7) 


Equations (3) through (7) are the governing relations used in this 
investigation of the SIP and ADIP numerical procedures. 


1 This assumption of 3T/3r = 0, r = 0, for all 4> values will be 
discussed in a subsequent section of this report. 
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jDUCIBtL^ OF THE 
■.tt pact 10 POOR 


1 1 1 1 1 F 


When the sphere is solid rather than hollow, a singularity exists at 
r ■ 0. At r • 0, the terms 


2k ar 
r 3r 


and 


r‘ sin <i 


3_ 

3$ 



are Indeterminate . ThrBe terms can be evaluated by using L'Hospltal'a rule 
(ref. 9). 2 


and 
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Likewise, a singularity exists at $ a 0, n in the term 


r‘ sin <> 


k ( Bln 


3T\ 
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"For the limits to exist, it is required that 
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Application of L'Hospltal'B rule to thiB term yield* 


r k 3 

/ u i_ * 22\1 
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Therefore, at the singularities, 
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The singularity at 


•an also be eliminated by approximating the center 


in a Cartesian form Tation. ’ This approach is discussed by Smith (ref. 10). 

A related question concerning a singularity for the cylindrical problcu is 
discussed by Albasiny (ref. 6 ). A third method of eliminating the singularity 
at r * 0 is to simply aBBume that a small but finite radius exists at the 
center; i.e., hollow-sphere approximation. All three approaches were examined 
in this investigation and will be discussed in a subsequent section of this 
report. 


J The singularities at r = 0 and $ = 0, it can also be eliminated by 
not specifying node points on the boundaries. 


6 


lMl’LICl DIFFERENCE EQUATION FORMULATION 


Equation (j) describes the heat transfer within a spherical -egion, and a 
solution Is achieved by approximating the partial derivatives with the use of 
sultub'e finite-difference expressions Involving the independent and dependent 
variables. The two-dimensional spherical region with an r, $ grid system 
imposed is shown in figure 2. An implicit central-difference equation for 
each grid point (i, j) within the specified region can be written as 



♦q ,M (9) 

where T* is the unknown temperature. Equation (9) can be written as 



Rewriting equation (10) yields 
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AT' + B T* + C T' + D T' + E T* = 

i,J i,J-l B i.J i-l«J i,J l+l ,J a i,ri,J+l 


\ ‘ 
x »c 


( 11 ) 


Equation (ll) has five unknown temperatures per grid point (l t j). The 
values of A, B, C t D, E, and Q are known on the basis of the previous 
time level and/or the previous iteration. A set of equations similar to 
equation (11 ) can be written for all i 4 j grid points within the region and 


ll 

on the boundaries. This matrix of equations can then be inverted to yield 
the unknowns , T* . For large systems of equations , this matrix solution can 
i »<J 

become very time consuming. 


NUMERICAL SOLUTION TECHNIQUE 


Stone (ref. 2) developed the SIP, an iterative method for solving sets 
of algebraic equations that occur for multidimensional systems. This method 
has been used successfully in solving problems involving two-dimensional, 
steady-state heat conduction, as well as multidimensional flow in a petroleum 
reservoir (ref. 3). The foundation of the SIP calculation method is based on 
the approximate factoring of the five-diagonal matrix (five nonzero elements 
in each row of matrix) generated by equation (ll) into three-diagonal upper 
and lower triangular matrices. The detailed mathematical reduction process 
required to derive the upper and lower triangular matrices is presented by 
Stone (ref. 2). The equations used in the SIP algorithm to solve for the 
unknown variable TJ , together with the boundary condition restrictions, can 

9 O 

be found in references 2 and 3. 


A second method used in the solution of equation (ll) Is the ADIP 
(Peaceman and Rachford (ref. 8)), which reduces the number of unknowns to 
three, as obtained for simple, one-dimensional problems. Basically, the ADIP 
solves the equations in one direction, with the dependent variable in the 
second dimension assumed constant over the time interval. As an example, 

consider equation (ll) for the first time step, in the $ direction."* 




CT! , 


DT' 

i+l.J 


= Q i,J ' AT i J-l “ ^'i.J+l 


( 12 ) 


The SIP boundary condition restrictions are illustrated in the appendix. 

Although not specifically pointed out, each time step is split into two 
parts. The first one-half time step is differenced implicitly in $ and 
explicitly in r, whereas the second one-half time step is differenced 
Implicitly in r and explicitly in $. 
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The temperatures 


i.J-1 


and T 




are known from the previous time 


step. Application of equation (12) to a grid network yields a tridiagonul 
matrix of unknown temperatures. The advantages of solving a tridiagonal matrix 
rather than a pentadi agonal matrix as generated by equation (ll) are evident. 


In addition to the previous two methods, the weighted average approach of 
Crank and Nicolson (ref. 11 ) waB used In conjunction with the SIP algorithm. 
The Crank-Nicolson (CN) modification is illustrated in the following 
application to equation (9). 
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(13) 
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Equation (13) can be rewritten into the fora of equation (11 ), which ia 
consistent with the f>IP formulation. For the CN method, 0 ia set equal to 

0,5 in equation (13).^ TMb method is designated simply as SIP/CN. 

Because the .'IP is an algorithm for solving a certain type of matrix, a 
table of geometrically arranged iteration parameters is normally erap'oyed to 
speed convergence. Weinstein et al. (ref. 3) recommended a geometrical itera- 
tion parameter defined by the relation 

JL 

1 •* Y ■ fl-Y \ ”“ 1 , m - 0, 1, ... (n - 1) (1U) 

m V max .1 

where y is the iteration parameter and n is the number of parameters 
(normally 1* to 10) in a cycle. The value of the iteration parameter lies 
between 0 and 1. For a neat-conduction problem with constant properties 
(ref. 2 ), 


v 

max. 


1 


min. 


?LArf . 

k.(Ar) 2 

. !_ 

k (r A$ ) 2 
r 


2(r A») 2 

k ( r A*) 2 

1 + “ — 

k^(Ar) 2 


(15) 


For this study, a maximum y value of 0.95 was used. A discussion of the 
physical and mathematical significances of the iteration parameter can be 
found in references 2 and 3. 


COMPARATIVE NUMERICAL RESULTS 


To study the effect of various boundary conditions on the relative 
accuracy of the solution techniques, the following three examples are 
7 

Co., idered. 


°Steen and Ali (ref. 5) used weighting values of 0.5 and 0.75. 

7 

Any consistent set of units can be used in these examples. In this 
study, absolute numerical values are used instead of dimensionless quantities, 
for comparison purposes. Numerical values of k = 0.8, C = 0 .I 4 , p * 130, 
and R = 1 have been used in these examples. 
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r 


Ca3e 1 - A homogeneous, two-dimensional, uniform surface temperature is 
specified. The boundary and initial conditions are 


T(R, ♦, t) * constant surface temperature 
|“(0, ♦, t) « 0 


uT 

gj(r, ♦, t) ■ 0, ♦ ■ 0, n 

q" • • 0 

T(r, ♦, 0) » T i 

This first example is for a sphere with a specified surface 

temperature. The surface boundary condition is such that 

T(R, #) * constant. An analytic solution is available (ref. 12). 

Case 2 - A homogeneous, two-demens ion, uniform heat generation is 
specified. 


T(R, +, t) - T(r, $, 0)) 


f“(0, ♦, t) * 0 

yn 

jrir, $, t) « 0, ♦ ■ 0, it 

q' "(r, ♦ , t) * constant 

Example 2 considers a sphere with uniform internal heat gener- 
ation. Both the initial and surface temperatures are set equal 
to zero. An analytic solution for this case can be found in 
reference 13. 

Case 3 - A homogeneous, two-dimensional, nonuniform surface temperature 
is specified. 
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T( Ft, ♦ , t) - r‘ 5 


|p<0, ♦ , t) • 0 

jy(r, *, t) - 0, ♦ ■ 0, w 

q' • » ■ 0 

T(r, *, 0) - 0 

This third example considers a sphere of unit radius R ■ 1, with 
the surface temperature specified as a function of cos $. An 
analytic solution for this case is presented in reference 7. 

It should be noted that cases 1 and 2 are one-dimensional problems; 
however, the numerical computations were performed with the use of a two- 
dimensional model. Case 3 is used to represent the accuracy of the numerical 
techniques for a two-dimensional problem with a zero temperature along 
$ * 5**, 7356 degrees and $ * 125.261*14 degrees for all values of r. For 
these cases, the results are given in terms of the difference between the 
temperature obtained by the exact solution and that obtained form the various 
numerical solution techniques. These results are called the temperature 
errors, defined as ( . As a convenient reference for comparing 

the numerical data, table I summarizes the various conditions used to generate 
the results given in tables II through IV. For example, numerical time-step 
effects can be studied by reference to tables 11(a) and 11(b). Tables 11(a) 
through 11(c) present a comparison of the numerical results for locations 
within a sphere with a specified constant surface temperature condition of 
r/R = 0 and r/R * 0.5. The temperature history at the geometrical center, 
r/R = 0, is of special interest because a discontinuity in equation (3) occurs 
at this location. A comparison of the results at r = 0 indicates that the 
SIP with CN modification (SIP/CN) with a geometrically variable y and the 
ADIP are the most accurate. A maximum temperature error of 30.631 degrees 
(3.96 percent error) occurred at a time unit of 10 after the start of the 
transient. Although no <J> variation in the surface temperature was 
specified, where a 4 > variation in the temperature was calculated, the error 
range in the 4 direction is shown. 

The standard SIP methods (constant and variable y) had the greatest 
absolute errors. Also shown in tables 11(a) through 11(c) are the hollow- 
sphere and rectangular approximation solutions used at r * 0. Again, a large 
absolute error was found for these two approximate solutions. The effect of 
the time step is shown in tables 11(a) through 11(c): reduction of the time 

step to At = 0.1 resulted in a significant reduction in the absolute error 
for all methods investigated. The effect of location, r/R * 0.5, on error 
again shows the SIP/CN (variable v) and ADIP methods to be the most accurate. 
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The effect of node size on accuracy can be seen by comparing the results 
of tables 11(b) and 11(d) for r ■ 0, At « 1.0. For a reduction of Ar ■ 0.10 
to Ar ■ 0.0!), the error with use of the SIP (y ■ constant) increased from 
30.631 to 92.310 degrees for At * 1.0. A similar increase in the temperature 

g 

error for the SIP/CN (y ■ vai table) and the ADIP was experienced. 

However, for a node reduction of Ar ■ 0.1 to Ar * 0.05, at a time step 
of At ■ 0.10, the absolute error decreased for both the SIP/CN and the ADIP 
methods. The results in tables 11(a) through 11(e) clearly Indicate the effect 
of time step and node size on numerical accuracy for the SIP approach. 

Tables 111(a) and Ill(b) present the numerical results for a sphere with 
internal heat generation. Once again, the SIP/CN (variable y) and the ADIP 
methods are the most accurate. For thiB particular case, a reduction in the 
time step of At * 1.0 to At * 0.1 resulted in a greater accuracy, in 
general, for the five methods, except for the SIP/CN (variable y) and the ADIP 
methods. 

Tables IV(a) through IV(g) present the numerical results for a Bphere with 
the surface temperature specified as a cos <> function. Table IV(a) 1 b the 
analytic solution as outlined in reference 7. For case 3, both an absolute 
error defined by 


T » T — T 

error exact cal. 


and a relative error defined by 


T 


rel 


T - T 

exact cal. 

T 

exact 


were used to evaluate the numerical procedures. Tables IV(b), IV(d), and 
IV(f) show the steady-state absolute error for ADIP, SIP/CN (a *= variable), 
and SIP (a = variable) /hollow-sphere approximation, respectively. Tables 
IV(c), IV(e), and IV(g) show the relative error for the respective methods. 

As expected, the least error occurs for r values near R = 1 and for $ 
values greater than 55 degrees and 125 degrees. As r approaches zero, the 
relative error increases quite rapidly. This same effect is observed as $ 
approaches 55.7** degrees and 125.26 degrees, where temperature 1 b zero for r 
values. These errors are a result of the assumption (e.g., eq. 5) used in 
numerical procedures and illust-ate the sensitivity when the solution is zero. 
Similar results are shown by Kee in reference 7, where the restrictions of 
equations (8a) and (8b) were not employed. 


A similar result was also noted by Barakat and Clark (ref. 1*4 ). 


CONCLUDING PttiAHKS 


Two basic numerical solutions (the strongly implicit procedure (SIP) and 
the alternating-direction implicit procedure (ADIP)) to the diffusion equation 
in spherical coordinates have been presented. The validity and accuracy of 
these solutions are demonstrated by comparing the results obtained therby with 
those of analytical solutions. Previous studies have shown that both methods 
comi>are favorably for the diffusion equation in Cartesian coordinates. The 
standard SIP appears to be slightly less efficient than the ADIP for the solid 
spherical problem studied in this investigation. This decrease in efficiency 
may be a direct result of the requirement that the dependent variable be 
calculated for the center of the sphere, where a discontinuity in the governing 
equation occurs. However, the Crank-Nicolson modification of the SIP gave 
essentially the same results as the ADIP for the cases studied. 

In conclusion, it should be mentioned that the SIP algorithm has been 
shown to be far superior to the ADIP for simulation problems involving multi- 
phase flow in porous media. It has been possible to obtain converged 
solutions to coupled systems of partial differential equations with the SIP 
when both the AI IP and successive over-relaxation procedures have failed. It 
should also be pointed out that a comparison in which a constant-property 
rectangular region was used showed the ADIP to be superior to the SIP, but the 
SIP - as shown more efficient for other rectangular cases involving property 
anis>. ophy and/or irregular boundary conditions. 


Lyndon B. Johnson Space Center 

National Aeronautics and Space Administration 
Houston, Texas, January 2U, 1977 
986-15-31-01-72 
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TABLE I.- .UMMARY OF CONDITIONS STUDIED 


r/H 

At 

Ar 

Table no. 

Remarks 

Case I 

0 

0.1 

0.10 

11(a) 


0 

1.0 

.10 

11(b) 


.5 

1.0 

.10 

11(c) 


0 

1.0 

.05 

11(d) 


0 

.1 

.05 

11(e) 


Case 2 

0.5 

1.0 

0.10 

m(i) 


.5 

.1 

.10 

Ill(b) 


Case 3 

— 

— 

— 

IV(a) 

Analytical solution 

— 

1.0 

0.025 

IV(b), IV(c ) 

ADIP 

— 

1.0 

.025 

IV(d), IV(e) 

SIP/CN (a = variable) 


1.0 

.025 

IV(f), IV(g) 

SIP (a *» variable)/ 
hollow-sphere 
approximation 
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Method 


Temperature error (T - T . ) a , deg, at 

£XHCT, CM • 

a time unit and corresponding T a t of - 



5 

10 

20 

50 


597.833° 

773.586° 

918.735° 

959.566° 

SIP (y - Y mB¥ ) 
max. 


2.606 

2.533 

1.815 

1.797 

0.051* 

SIP (y “ variable) 

-l*.07l* 

2.1U6 

1.637 

.01*9 

SIP/CN (y - Y ) 

max . 

-2.978 

-3.006 

.563 

.5»*0 

.735 

.731 

.025 

SIP/CN (y * variable) 

-3.099 

.1*28 

.685 

.023 

ADIP 

-3.037 

.1*25 

.680 

.026 

SIP ( Y - Y m „ ¥ )/ 
max . 

hollow -sphere 
approximation 

2.973 

7.805 

1.969 

-1.381* 

SIP « Y_ )/ 

max. 

rectangular 

approximation 

3.862 

9.017 

3.355 

.071 


* -xact temperature; T . ■ calculated temperature. 

6X&CL Cfti * 
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Method 


SIP (y ■ y ) 

'max. 

SIP (y * variable) 

SIP/CN (y * Y ) 

'max. 

SIP/CN (y * variable) 


ADIP 

SIP ( Y - Y_„ )/ 
max • 

hollow-sphere 

approximation 

SIP (y = Y_ )/ 

max • 

rectangular 

approximation 


Temperature error (T - T , ) a , deg, at 

0X&C v* CU1 • 

a time 

unit and corresponding T of - 

5 

10 

20 

50 

597.833° 

773.586° 

918.735° 

959.566° 

-3.09*4 

30.631 



-.198 

25.886 


.5*47 

-10.293 

15.585 
15. *485 

10.309 

10.283 

.33*4 

.1**1 

9.986 

3.91*4 

.09 

*4.255 

*4.886 

2.959 

.07 

-.821 

.30*4 

.*460 

.*455 

.017 

1.095 

-1.13 

-.306 

-.002 

8.019 

35.133 

17.083 

-.875 

8.525 

36.09*4 

18. *4 56 

-.1*45 


















TABLE II.- Continued 


(c) r/B ■ 0.5i At ■ 1.0; Ar ■ 0.10 


Method 

Temperature error (T rjtact - T cal> ) a » de *» at 


a time 

unit and corresponding T t , X(lct * 


5 

10 

20 

50 


703.751° 

81*0.065° 

933.727° 

959.724° 

S Ip (Y - Y m „, ) 
max • 

15.22*7 

15.**72 

18.891 

18.816 

9.563 

9.1*86 

0.339 

.338 

SIP (y ■ varieble) 

10.301 

13.576 

6.583 

6.580 

.213 

SIP/CN (y ■ Y ) 

max • 

1.150 

1 . 112 * 

2.189 

2.162 

1.298 

1.259 

.041 

SIP/CN (y - variable) 

-.880 

.165 

.289 

.011 

ADIP 

3.118 

-.867 

-.285 

-.002 

SIP ( Y “ W >/ 

max • 

hollow-sphere 

approximation 

15.22*7 
lb . 288 

18.891 

18.309 

9.563 

9.361* 

.338 

SIP )/ 

max . 

rectangular 

approximation 

15.22*7 

lJ *.288 

18.891 

18.816 

9.563 

9.361* 

.338 


R T ^ , » exact temperature; T cal * calculated temperature. 
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TABLE II.- Continued 


(d) r/R ■ 0; At ■ 1.0; Ar ■ 0.05 
Method Temperature error (T • T . )*, deg, at 

6XHC ti CM • 


a time unit and corresponding T Rf>t - 



5 

10 

20 

50 


597.833° 

773.586° 

918.735° 

959.566° 

SIP (y • y ) 

28.523 

33.555 

92.310 

80.253 

55.256 

50.11*3 

2.751 

2.579 

SIP (y ■ variable) 

2.065 

l*.6o6 

32.138 

28.219 

16.972 

15.838 

.536 

.517 

SIP/CN (y - Y ) 

max . 

18 .U 5 U 

30.1*51 

5»* • 256 
26.880 

25.711 

16.711* 

1.251 

.921* 

SIP/CN (y - variable) 

2.568 

6.652 

10.233 

1.126 

5.86y 

1.013 

.203 

.288 

AJJIP 

3.61*0 

-1.951 

-.869 

-.013 


^exact " exac * temperature; T cr1 ■ calculated temperature. 
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TABLK II.- Concluded 


(e) r/K ■ 0; At - 0.1; Ar ■ 0.05 


Method 

T*i«r«tu» -rror dm. «t 

a time unit nr»d corresponding t of - 


5 

10 

20 

50 


597.833° 

771.586° 

918.735° 

959.566° 

BIT (y - y___ ) 

max . 

l.i*6o 

1.572 

6.112 

6.206 

2.881* 

2.900 

0.073 

SIP (y ■ variable) 

-1.799 

-1.790 

1.088 

1.887 

1.123 

1.12U 

.028 

SIP/CN (y ■ Y ) 

'max. 

.152 

.1*17 

1.301* 

1.583 

.723 

.756 

.017 

SIP/CN (y * variable) 

-.768 

.112 

.168 

.003 

ADIP 

-.699 

.110 

.166 

.008 


*T 

exact 


exact temperature; T 

cal. 


calculated temperature. 
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tab>j: III.- CASE 2 


(a) r/R ■ 0.5; At * 1.0; Ar ■ 0.10 


Method 

Temperature error (T . - T , ) a , 

exact cal . ' 

deg, at a time unit and corresponding 

T of - 

exact 


10 

2.1*l81*9 J 

30 

3.0911° 

50 

3.12337° 

SIP W ) 

mux • 

0.12552 

.121*07 

0.02150 

.02117 

0.00206 

.00203 

SIP (y * variable) 

.08681 

.08676 

. 011*11 

.00127 

SIP/CN (y » Y mBY ) 
max . 

.02132 

.02030 


.00026 

SIP/CN (y s variable) 

.00679 

.00099 

.00008 

ADIP 

.00038 

.00010 

0 


®T = exact temperature; T , * calculated temperature, 

exact cal. 


reproducibility of the 
#RWCNAL PACS DB POOH 
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"ABLE III.- Concluded 


(b) r/H - 0.5; At ■ 0.1; Ar - 0.10 


Method 

Temperature error (T - T . ) ft . 

exact cal. * 

detf* at a time unit and correspond Inf: 

T . of - 

exact 


10 

?.l4l8U9° 

— 

30 

3.0911° 

50 

1. 12337° 

SIP (y « y ) 

max. 

0.01796 

.01788 

0.00266 

0.00021 

3IP (y * variable) 

.01689 

.00250 

.00020 

SIP/CN (y * y ) 

max. 

.00910 

.00908 

.0013*4 

.00011 

GIP/CN (y « variable) 

.0088 

.00129 

.00010 

ADIP 

.0087*4 

.00129 

.00011 


',. v „ * exact temperature; = calculated temperature. 








TABLE IV.- CASE 3 


(a) Analytic solution 


Angular 


Tenperature (T), 

deg, at a radial position ( 

r) of - 


position 
(♦), deg 

1.00 

0.75 

0.50 

Q.25 

0.15 

0.10 

0.05 

— 

0 

1.0000 

0.56250 

0.250000 

0 .062500 

0.02250 

0.01000 

0.002500 

16.1*6150 

.849582 

.477890 

.212395 

.0530989 

.0191156 

.00849582 

.00212396 

36.9231 

.458663 

.257998 

.114666 

.0286664 

.0103199 

.00458663 

.00114666 

55.38U6 

-.0159537 

-.00897393 

-.00398841 

-.000997103 

-.000358957 

-.0001595?" 

-.0000398841 

73.81*62 

-.383892 

-.215940 

-.0959731 

-.0239933 

-.00863758 

-.0038389: 

-.000959731 

83.0769 

-.478206 

-.268991 

-.119552 

-.0298879 

-.0107596 

-.00479206 

-.00119552 

96.9231 

-.478206 

-.26899: 

-.119552 

-.0298879 

-.0107596 

-.00478206 

-.00119557 

106.151*0 

-.383893 

-.21591*0 

-.0959731 

-.0239933 

-.00863758 

-.00383893 

-.000959732 

124.615 

-.0159538 

-.00897404 

-.00398846 

-.000997115 

-.000358961 

-.000159538 

-.0000398846 

143.077 

.458663 

.257998 

.114666 

.0286664 

.010319° 

.00458663 

.00114666 

161.538 

.849582 

.477890 

.212395 

.0530989 

.0191156 

.00849582 

.00212395 

180.00 

1.0000 

.56250 

.250000 

.062500 

.02250 

.01000 

.002500 


ru 

V/! 



ro 

CN 


TABLE IV.- Continued 


(b) AD IP - absolute error at steady-state 


^ Ingular 

Temperature absolute error (T - T ) a , 

6xsc L cai • 

deg, at a radial position ( 

r) of - 

(♦). deg 

1.00 

0.75 

0.50 

0.25 

0.15 

0.10 

0.050 

0 

0 

-0.000708 

-0.00086U 

-0.0006358 

-O.OOOL991 

-0.0001*358 

-0.000391*90 

18.U6150 

0 

-.000632 

-.000779 

-.0005919 

-.0001*773 

-.0001*21*02 

-.00038105 

36.92 31 

0 

-.0001*11* 

-.00051*1* 

-.0001*71*9 

-.OOOU196 

-.00039289 

-.00037093 

55.3846 

0 

-.00011*885 

-.00026036 

-.000332609 

-.00031*91*51* 

-.00035506 

-.000358630 

73.81*62 

0 

.000057 

-.0000398 

-.0002223 

-.00029508 

-.00032572 

-.00031*9093 

83.0769 

0 

.00011 

.000016 

-.000191* 

-.0002811 

-.0003182 

-.00031*6652 

96.9231 

0 

.00011 

.000016 

-.000191* 

-.0002811 

-.0003192 

-.00031*6652 

106.151*0 

0 

.000057 

-.0000398 

-.0002223 

-.00029507 

-.00032573 

-.00031*9091* 

12l*.6l50 

0 

-.00011*88 

-.00026031* 

-.000332612 

-.00031*91*56 

-.00035507 

-.000358633 

11*3.077 

0 

-.0001*11* 

-.00051*1* 

-.0001*71*9 

-.OOOU196 

-.00039290 

-.00037093 

161.538 

0 

-.000632 

-.000779 

-.0005919 

-.0001*77'' 

-.000171*03 

-.00038107 

180.00 

0 

-.000708 

-.000861* 

-.0006359 

-.0001*991 

-.0001*358 

-.000381*91 


exact 


T 


cal 


= exact temperature; 


calculated temperature. 


TABL£ IV.- Continued 


(c) ADIP - relative error at steady-state 


Angular 
position 
(♦), deg 

Temperature relative error ((T . - T , )/T l a at f 

exact cal . exact 

1 radial position (r) of - 

1.00 

0.75 

0.50 

0.25 

0.15 

0.10 

0.050 

0 

0 

-0.0012587 

-0.0031*56 

-0.010173 

-0.022182 

-0.01*3580 

-0.15396 

18.1*6150 

0 

-.0013225 

-.003667 

-.01111*7 

-.021*969 

-.01*9910 

-.1791*10 

36.9231 

0 

-.001601*7 

-.001*71*1*5 

-.016568 

-.01*0659 

-.085660 

-.3231*9 

55.38U6 

0 

.016587 

.065279 

.33358 

.97 353 

2.2256 

8.9918 

73.8U62 

0 

-.00026396 

.0001*11*70 

.0092651 

.031*162 

.081*81*7 

.36371* 

83.0769 

0 

-.0001*0891* 

-.00013383 

.0061*909 

.026126 

.06651*0 

.28996 

96.9231 

0 

-.0001*0891* 

-.00013383 

.0061*909 

.026126 

.06651*0 

.28996 

106.151*0 

0 

-.00026396 

.0001*11*70 

.0092651 

.031*161 

.061*81*- 

.36 371* 

12U.6150 

0 

.016581 

.065273 

.33357 

• 97352 

2.2256 

8.9916 

11*3.077 

0 

-.001601*7 

-.001*71*1*5 

-.016568 

-.01*0659 

-.085662 

-.3231*9 

161.538 

0 

-.0013225 

-.0036677 

-.01111*7 

-.021*969 

-.01*9910 

-.1791*2 

180.00 

0 

-.0012587 

-.0031*56 

-.010171* 

-.022182 

-.01*3580 

-.15396 


®T 

exact 


T 


cal . 


= exact temperature ; 


= calculated temperature. 



ro 

CD 


TABLE IV.- Continued 


(d) SIP/CN (a = variable) - absolute error at steady-state 


Angular 

Temperature absolute error ~ 

deg, at a radial position ( 

r) of - 

( $ ) , deg 

1.00 

-.75 

0.50 

0.25 

0.15 

0.10 

0.050 

0 

0 

-0.000716 

-0.000865 

-O.OOO6096 

-0.00011*65 

-0.0003796 

-0.00033051 

18.U6150 

0 

-.00061*0 

-.000779 

-.0005912 

-.0001793 

-.00012905 

-.00036815 

36.9231 

0 

-.0001*18 

-.00051*1* 

-.0001*71*7 

-.0001*187 

-.00039283 

-.00037706 

55.38U6 

0 

-.00011*826 

-.00026011 

-.000332297 

-.00031*9282 

-.000353759 

-.00036516 

73.8162 

0 

.00006 

-.0000395 

-.000222 

-.00029188 

-.00032509 

-.00035H08 

83.0769 

0 

.000115 

.000017 

-.0001937 

-.0002808 

-.00031788 

-.0003167 

96.9231 

0 

.000115 

.000017 

-.0001937 

-.0002608 

-.00031826 

-.000313671 

106.151*0 

0 

.00006 

-.OOOC395 

-.000222 

-.000291*72 

-.0003260 

-.00031117 

12U.615 

0 

-.00011*832 

-.0002602 

-. 000332*23 

-.00031*9281 

-.00035566 

-.0003517636 

113.077 

0 

-.0001*18 

-.00051*1* 

-.0001*71*7 

-.0001198 

-.00039215 

-.0003615 

161.538 

0 

-.000639 

-.000779 

-.0005922 

-.0001752 

-.000H67? 

-.00039205 

180.00 

0 

-.000716 

-.000865 

-.0006958 

-.0006199 

-.0005671 

-.00051088 


®t „ 

exact 


= exact temperature ; 


T 1 = calculated temperature. 

C&l • 



TABLE IV.- Continued 


(e) SIP/CN (o = variable) - relative error at steady-state 


Angular 

Temperature relative error ( (T gxact - T cal . )/T exa ct )& at a 

radial position (r) of - 

(♦), deg 

1.00 

0.75 

0.50 

0.25 

0.15 

0.10 

0.050 

0 

0 

-0 .0012729 

-0.00346 

-0.0097536 

-0 .019844 

-0.03796 

-0.132216 

18.U6150 

0 

-.0013392 

-.0036677 

-.011134 

-.025074 

-.05050 

-.17333 

36.9231 

0 

-.0016202 

-.0047442 

-.016559 

-.040572 

-.085647 

-.32883 

55.3846 

0 

.016521 

.065216 

.33326 

.97305 

2.2174 

9.15553 

73.8462 

0 

-.00027785 

-.00041157 

.0092526 

.034139 

.084682 

.36615 

83.0769 

0 

.00042752 

-.0001422 

.0064809 

.026098 

.066473 

.28999 

96.9231 

0 

.00042752 

-.0001422 

.0064809 

.026098 

.066552 

.28747 

106.1540 

0 

-.00027785 

-.00041157 

.0092526 

.034121 

.084919 

.35892 

124.615 

0 

.016528 

.065236 

.33338 

.97304 

2.2293 

8.81953 

143.077 

0 

-.0016202 

-.0047442 

-.016559 

-.040679 

-.085564 

-.31814 

161.538 

0 

-.0013371 

-.0036677 

-.011153 

-.024859 

-.049050 

-.18459 

180.00 

0 

-.0012729 

-.00346 

-.011133 

-.027551 

-.0567U 

-.204352 


= exact temperature; T , = calculated temperature, 

exact cal . 




uo 
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TABLE IV.- Continued 


( f ) SIP/hollow-sphere approximation (a = variable) - absolute error at steady-state 


Angular 

Temperature absolute error - T pal ) a . 

deg, at a radial position ( 

r) of - 

1 L 1UU 

(♦) , deg 

1.00 

0.75 

0.50 

0.25 

0.15 

0.10 

. 

0 

0 

- 0.000706 

-O.OCO 862 

- 0.0006318 

-O.OOOU 9 U 2 

-0.0001*307 

-0.0003791*5 

18 .U 6150 

0 

-.000631 

-.000777 

-.0005878 

-.OOOU 72 I* 

-.0001*1^36 

-.00037561 

36.9231 

0 

-.0001*13 

-.00051*2 

-.0001*707 

-.0001*11*7 

-.00038771* 

-.00036551* 

55.38U6 

0 

-.00011*737 

-.00025767 

-.00032823 

-.0003U577 

-.00031.9921* 

-.0003532921 

73.8462 

0 

.000058 

-.0000369 

-.0002179 

-.00029017 

-.0003206 

-.00031*3791* 

83.0769 

0 

.000112 

-.000019 

-.0001896 

-.0002762 

-.00031307 

-.00031*1361* 

96.9231 

0 

.000112 

-.000019 

-.0001896 

-.0002762 

-.00031307 

-.00031*1367 

106.151*0 

0 

.000058 

-.0000369 

-.0002179 

-.00029017 

-.0003206 

-.000 31*3803 

12 U .6150 

0 

-.0001U733 

-.00025766 

-.000328325 

-.0003l»l*587 

-.00031*9936 

-.0003533106 

11*3.077 

0 

-.0001*13 

-.00051*2 

-.0001*707 

-. 000 U 1 L 8 

-.00038777 

-.00036557 

161.538 

0 

-.000631 

-.000777 

-.0005878 

-.0001*725 

-.0001*189 

-.00037567 

180.00 

0 

-.000706 

-.000862 

-.0006318 

-.0001*91*3 

-.0001*307 

-.0003795 


exact 


cal , 


= exact temperature ; T 


= calculated temperature. 
















TABLE IV.- Concluded 


(g) SIP/hol low-sphere approximation (a = variable) - relative error at steady-state 


Angular 
position 
(4i), deg 

Temperature relative error ( ( T exact “ T ca l . ^ ^exact ^ &at a 

radial position (r) of - 

1.00 

0.75 

0.50 

0.25 

0.15 

0.10 

0.050 

0 

0 

-0.0012551 

-0.003UU8 

-0.010109 

-0.021961* 

-0.04307 

-0.15178 

18.6U15 

0 

-.001320U 

-.0036583 

-.011070 

-.021*713 

-.01*9302 

-.17681* 

36.9231 

0 

-.0016008 

-.001*7268 

-.0161*20 

-.01*0181* 

-.091*537 

-. 31879 

55.381*6 

0 

.0161*22 

.061*605 

.32918 

.95991* 

2.1931* 

8.8580 

73.8U62 

0 

-.00026859 

.00038UU8 

.0090817 

.039591* 

.083513 

.35822 

83.0769 

0 

-.0001*1637 

.00015893 

.00631*37 

.025670 

.0651*68 

.28551* 

96.9231 

0 

-.0001*1637 

.00015893 

.00631*37 

.025670 

.0651*68 

.28551* 

106.15U0 

0 

-.00026859 

.00038UU8 

.0090817 

.033591* 

.083513 

.35823 

12U.6150 

0 

.0161*17 

.061*601 

.32927 

.95996 

2.1931* 

8.8583 

11*3.077 

0 

-.0016008 

-.001*7268 

-.0161*20 

-.01*0191* 

-.O8U5U 

-.31881 

161.538 

0 

-.0013201* 

-.0036583 

-.011070 

-.021*718 

-.01*9302 

-.17687 

180.00 

0 

-.0012551 

-.0031*1*8 

-.010109 

-.021969 

-.01*307 

-.15178 


^ = exact temperature; 

exact 


T = calculated temperature. 

C&l • 


UJ 




/ 



33 



APPENDIX 


BOUNDARY CONDITION RELATIONS WITH USE OF 
THE STRONGLY IMPLICIT TECHNIQUE 


The transient-heat-conduction equation in spherical coordinates 
(eq. (3)) is put in finite-difference form at r = 0 as an illustration of 
the SIP boundary condition requirements. 


SOLID SPHERE 


The SIP boundary restrictions require that 

A i,o’ E i,R * 0 for * ' °- ” 

Vo,j'W“ for r = 0,R 


(16) 

(17) 


As an illustration of these boundary conditions, consider the singularity 
located at the geometrical center, r = 0, of the sphere. 

By employing the boundary conditions represented by equations (8a) and 
(8b), equation (3) becomes 


DC 21 = ^T 

pc P at ^ 9r 2 + * 


(18) 


which can be written in finite-difference form as 


, - 6T! , + 3V . /T ' , - T \ 

..1 \ i,.1 + l 1*J~1/ + a « » * = pC ( ^ i-l-i j 

(fir ) 2 p \ At / 


(19) 


Employing equation (5) yields 


m — m 

i,J+l i.J-l 


( 20 ) 


Then, equation (19) can be written as 
PC 

2.1 tpt 

» j 


[ 6k. . pC "1 
LuL + I m i 

(Ar) 2 At J 1 ’ ; 




(Ar) 


2 i ,J+1 


= -q 


f » t _ 


At A i,j 


3H 


( 21 ) 


Comparing equation (21) with equation (11 ) yields 


where 


and 




♦ E 


i ,J T i .>1 




(22) 


(23a) 


(23b) 

(23c) 


(23d) 


HOLLOW-SPHERE APPROXIMATION 


This approximation assumes that a small but finite radius (r Q ) can be 

used to represent the geometrical center. To illustrate this boundary 
condition, consider the location 

r * r 0 (f 0 3 0.01 Ar ) * ♦ * 0 

By employing the boundary conditions represented by equation (8c) and 


equation (3) becomes 
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With the aastimptlons that at 
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equal on (2M can be formulated in terms of equation (il) as 
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